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Abstract 



The mass spectrum of S and P-wave mesons containing a single heavy quark 
is computed in the quenched approximation, using NRQCD up to third or- 
der in the inverse heavy quark mass expansion. Previous results found third 
order contributions which are as large in magnitude as the total second order 
contribution for the charmed S-wave spin splitting. The present work consid- 
ers variations such as anisotropic lattices, Landau link tadpole improvement, 
and a highly-improved light quark action, and finds that the second order 
correction to the charmed S-wave spin splitting is about 20% of the leading 
order contribution, while the third order correction is about 20%(10%) for 
D* — D(D* — D s ). Nonleading corrections are very small for the bottom 
meson spectrum, and are statistically insignificant for the P-wave charmed 
masses. The relative orderings among P-wave charmed and bottom mesons, 
and the sizes of the mass splittings, are discussed in light of experimental data 
and existing calculations. 
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I. INTRODUCTION 



Calculations in QCD can be performed numerically on a discrete space-time lattice, 
provided the lattice spacing is small enough to accommodate all of the revelant physical 
distance scales. In the presence of a heavy quark the lattice spacing must be small relative to 
the inverse quark mass, resulting in large computational requirements, unless an appropriate 
"effective theory" is used. In particular, for the case of a hadron containing exactly one 
heavy quark, the dynamics can be expanded in powers of the inverse heavy quark mass 
using the well-established technique of heavy quark effective theory [|IJ or nonrelativistic 
QCD (NRQCD) With the heavy quark expansion, the lattice spacing can be much 

coarser and the computational requirements are correspondingly smaller. 

In the present work, two issues will be addressed within quenched lattice NRQCD. Firstly, 
is the charm quark sufficiently heavy to permit the use of lattice NRQCD for charmed meson 
spectroscopy? Secondly, what are the mass splittings (magnitude and sign) between P-wave 
mesons containing a single heavy quark? 

The first issue is clearly of interest due to the computational efficiency of lattice NRQCD. 
The heavy quark expansion is known to work well for bottom mesons [fH|, but previous 
research has demonstrated that 0(1/M 2 ) and 0(1/M 3 ) contributions can be comparable in 
magnitude for the S-wave spin splitting of charmed mesons. [|| The present work provides an 
extension of this investigation by considering new simulations that incorporate a number of 
changes in method. For example, Ref. || used the fourth root of an elementary plaquette to 
define the tadpole improvement factor, Uq, whereas the present work uses the mean link in 
Landau gauge. Since the Landau definition is advantageous in other contexts ||, including 
the velocity expansion for the charmonium spectrum of lattice NRQCD 0, it might be 
expected to improve the convergence of the heavy quark expansion for the charmed heavy- 
light spectrum as well. Also, the light quark was described by the Sheikholeslami-Wohlert 
action in Ref. f4|, but the present work uses a D234 action || which has smaller lattice 
spacing errors classically. The simulations in the present work differ from those of Ref. 
Q in various other ways as well, including Dirichlet versus periodic boundary conditions 
for light quark propagation, differing discretizations for heavy quark propagation, and the 
introduction of an anisotropic lattice with a smaller temporal lattice spacing than spatial 
spacing. Despite these modifications, the present simulations produce a conclusion similar to 
that of Ref . |4j]: the 0(1/M 2 ) and 0(1/M 3 ) contributions are comparable in magnitude for 
the S-wave charmed spin splitting; in the present work, each is roughly 20% of the 0(1/M) 
result for D* — D, though the 0(1/M 3 ) contribution is closer to 10% for D* — D s . 

The second issue under discussion relates to the spectrum of P-wave mesons containing a 
single heavy quark. The relative orderings of the P-wave bottom mesons have only recently 
come under direct experimental scrutiny [pHT^], and the complete picture is not yet clear. 
Meanwhile theoretical predictions differ from one another even at a qualitative level. The 
traditional expectation of a hydrogen-like spectrum that arises from a number of model 
calculations |fH|-|T%[ has been questioned long ago by Schnitzer [TBI and very recently by 



Isgur |^D[ and by Ebert, Galkin and Faustov plf . 

The calculation is difficult within lattice QCD because the P-wave splittings are not 
large in comparison to the typical scale of nonperturbative QCD, and because of operator 
mixing for the pair of P-wave states having J = 1. Some previous attempts have been made 
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26|. Unfortunately, the uncertainties are often substantial, and results are not always 
as consistent with one another as might have been hoped. The present work represents a 
further comment on this situation. In particular, the D* 2 — D* mass splitting is found to 
be positive as in the traditional hydrogen-like ordering, and to be substantially less than 
100 MeV, while the D\ — D$, B* s2 — B* Q and B\ — Bq splittings are even smaller. These 
splittings are consistent with a number of model calculations, but are somewhat smaller 
than the lattice NRQCD calculation of Ref. |?3 . 



II. ACTION 

The lattice action has three terms: gauge action, light quark action and heavy quark 
action. The entire action is classically and tadpole-improved with the tadpole factors, u s 
and u t , defined as the mean links in Landau gauge in a spatial and temporal direction, 
respectively. 

The gauge field action is 



Sa(U) = f 



-E 



ps 



u i u t pt 



-ReTrC/, 



1 



1. 



20uM rts 



-ReTrC/, 



pt 



20u% 



20uju 2 t rst 



Ei 



1 



-ReTrf/ r t s 



Ei 



-ReTrC/ re 
3 



-ReTrC/,.. 



(1) 



where £ = a s /a t . The subscripts "ps" and "rs" denote spatial plaquettes and spatial pla- 
nar 1x2 rectangles respectively. Plaquettes in the temporal-spatial planes are denoted by 
"pt", while rectangles with the long side in a spatial (temporal) direction are labeled by 
"rst" ( "rts" ) . The leading classical errors of this action are quartic in lattice spacing. 

For light quarks, a D234 action || is used with parameters set to their tadpole-improved 
classical values. Its leading classical errors are cubic in lattice spacing. 
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D 2t (x) = ij){x){\ - ^ 4 )U 4 (x)U 4 (x + i)tjj(x + 2t) 
+ $(x + 2t)(l + i±)U\(x + t)Ul(x)il>(x). 



1 

2i 



-Im(Tr^(x)) 



(6) 
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+Ul(x - u)U^x - v)U v (x + fx- v)UUx) . 



(8) 



The heavy quark action is NRQCD 
function propagation: 



], which is discretized to give the following Green's 
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with n = 5 chosen for this work. Separation of the Hamiltonian into two terms, H = 
Ha + Hb, is important for ensuring stability of the discretization. For example, recall the 
discussion in Ref. of a large nonzero vacuum expectation value for the term containing 
Cio in the Hamiltonian (see Eq. (|l^)). This issue will be discussed further in Sec. [IV]. 

The following Hamiltonian, written in terms of the bare heavy quark mass M, is complete 
to 0(1/M 3 ) in the classical continuum limit [ 27| : 
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The coefficients of the Hamiltonian are chosen so the dimensionless parameters, c$, are 
unity at the classical level. As will be discussed below, computations have been performed 
with the Cj set to unity or zero in various combinations, including separate computations at 
0(1/M), 0(1/M 2 ) and 0(1/M 3 ) to allow discussions of convergence for the 1/M expansion. 
Throughout this work, H is always placed in Ha of Eq. @ and all of the remaining terms 
except the c\q term are only placed in Hb- The difference between having the cio term in 
Ha or Hb will be discussed explicitly, since it has the nonzero vacuum expectation value. 



4 



A tilde on any quantity indicates that the leading discretization errors have been removed. 
In particular, 

Ei = F 4l , (16) 
B% = -Eij k Fj k , (17) 

where || 

5 1 1 

F^(x) = -F^(x) - —U^(x)F^(x + jj)Ul(x) - Q-^Ufa - jl)F^{x - fi)U^(x - p) 

-(/i - u). (18) 
The various spatial lattice derivatives are defined as follows: 

a s A,G(x) = —[Ui(x)G(x + i) - Uj{x - t)G(x - i)}, (19) 

a s A^G(x) = ^^-G{x + z) - G(x), (20) 

a s A { r ] G{x) = G(x) - U ^ X ~^ G{x - z), (21) 

u s 

a 2 A? ] G(x) = ^-G{x + z) - 2G{x) + U ^ X ~^ G(x - z), (22) 

A, = A, - IaW^A^, (23) 
A( 2 ) = EAS 2) , (24) 

i 

A( 2 ) = A( 2 '-^AW, (25) 
12 V ; 

,(2)^ 2 



A (4) = EK } ) • (26) 

i 

This NRQCD action has quadratic classical lattice spacing errors. 



III. METHOD 



The data presented here come from 2000 gauge field configurations on 10 x30 lattices 
at /3 — 2.1 with a bare aspect ratio of £ = a s /a t = 2. Two light quark masses are used, 
corresponding to n — 0.23 and 0.24. Fixed time boundaries are used for the light quark 
propagators so they fit naturally into a heavy-light meson, since the NRQCD heavy quark 
propagator is also not periodic in the temporal direction. 

A calculation of the string tension from these gauge field configurations provides a de- 
termination of the renormalized anisotropy: 

f = a s /a t = 1.96(2). (27) 
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TABLE I. Heavy-light meson creation operators. 



2S+1 Lj 0(f) 

3 Si (0,a<) 

x Pi (OA) 

3 Pi (o Ao-j - A i^) 

3 P 2 (0,Vi - Aj-Oj-) or 



Using light quarks only, the lightest pseudoscalar and vector meson masses are easily 
obtained from local creation operators. By linear interpolation and extrapolation in 1/k, 
the critical (ac c ) and strange (k s ) hopping parameters and the temporal lattice spacing are 
found to be 

k c = 0.243025(41), (28) 

_ /0.2344(H) from m , 
s \0.2356(3) fromm K , 1 > 

at = 0.1075(23)fmfromm p . (30) 

For the hopping parameters used in explicit computations, k = 0.23 and 0.24, the ratio 
of pseudoscalar to vector meson masses is "m ff /m p " = 0.815(3) and 0.517(8) respectively. 
No exceptional configurations were encountered at these k values. One might expect a 
systematic uncertainty on a t to account for deviations from the linear relationship between 
a t m p and 1/k. The uncertainty is presumably a few percent, but cannot be estimated using 
only the two k values studied here. 

A heavy-light meson is created by the following operator, 

]TQt (fMf)r(f)g(f ) 5 ( 31 ) 

X 

where Q(x) is given in Table | and the smearing operator is 

T(x) = [l + c s A^(x)} n \ (32) 

All plots shown here use (c s ,n s ) = (0.15, 10) at the source and a local sink. The source is 
fixed at timestep 4, which is a distance 3a t from the lattice boundary. 

Because NRQCD is an expansion in the inverse bare heavy quark mass, all meson mass 
differences can be obtained from correlation functions at p = 0, but the absolute meson 
mass itself cannot be obtained directly. One way to determine the mass is to compute the 
change in energy when a meson is boosted, 

P 2 

E V -E G = - J — . (33) 
P 2M kin 1 ; 

This defines the kinetic mass, M^ in , which is interpreted as the meson's physical mass. For 
the present work, E p is computed only for the 1 Sq state, with p = (0, 0,2n/L s ) where 
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FIG. 1. Effective mass plots for the ground state heavy-light meson at rest with k = 0.23, and 
terms up to 0(l/(a s M) 2 ) retained in the Hamiltonian. Only Ho is in Ha of Eq. (|9|). Solid circles, 
solid squares and open circles are for a s M = 1.2, 1.5 and 5.0 respectively. Statistical errors are 
smaller than the plotted symbols. 



N s a s is the spatial extent of the lattice in physical units. Solving for the kinetic mass 



gives 



Mi 



2tt 2 



kin 



N^at[ot(E p -E )) 



(34) 



Some justification for Eq. (|33"D comes from consideration of the next correction term, giving 



E p — E 



2M, 



kin 



8M k y 



(35) 



We have verified that the extra term shifts the kinetic mass by an amount which is smaller 
than the uncertainties for every value of Mk in reported in this work. 

In the case of S-wave mesons, a plateau containing ample timesteps is clearly evident in 
all effective mass plots; some examples are shown in Fig. p]. A more detailed discussion of 
the P-wave plateaus is deferred to Sec. |V|. In this paper, the plateau region for each mass 
is defined by the maximum value of 
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TABLE II. Examples of plateaus defined by maximization of Q. Notice that they all have the 
desired feature that Q > 0.1. Insensitivity of the fit parameters to the precise plateau boundaries 
has been verified. 





So 


%(p) 


3 Si 


'Pi 


3 Po 


3 A 


3 P 2 


a s M = 1.5, k 

Q 


= 0.23 
(15,22) 
0.45 


(16,22) 
0.13 


[17,22) 
0.27 


(13,22) 
0.60 


/in O 1 ^ 

(12,22) 
0.21 


( 1 /I 00\ 

(14,22) 
0.97 


(12,22) 
0.92 


a s M = 1.5, k 

Q 


= 0.24 
(15,20) 
0.22 


(16,20) 
0.14 


(15,20) 
0.61 


(12,20) 
0.69 


(11,20) 
0.26 


(14,20) 
0.86 


(13,20) 
0.43 


a s M = 6, k = 

(^mirij ^max) 

Q 


0.23 
(17,22) 
0.88 


(16,22) 
0.29 


(16,22) 
0.67 


(13,22) 
0.15 


(12,22) 
0.69 


(13,22) 
0.65 


(12,22) 
0.50 


a s M = 6, k = 

(^miri) ^max) 

Q 


0.24 
(16,20) 
0.51 


(14,20) 
0.16 


(16,20) 
0.77 


(13,20) 
0.35 


(12,20) 
0.10 


(13,20) 
0.12 


(12,20) 
0.93 






Q = 


T(N/2 - 
T(AT/2 


l,X 2 /2) 
-1,0) 






(36) 


where 




T(a, x) ■ 


= / dtt' 


1-1 exp(-t), 






(37) 



and N is the number of timesteps in the proposed plateau region, \ 2 is obtained from a 
single exponential fit to each meson correlation function. A correlated fit is done with the 
covariance matrix inverted by singular value decomposition. Statistical uncertainties are 
obtained from the analysis of 5000 bootstrap ensembles. All plateaus are ended at timestep 
22(20) for simulations with k = 0.23(0.24), except for an infinitely-heavy quark due to 
excessive noise at these larger times. The examples in Table demonstrate the quality of 
the fits. 

IV. S-WAVE SPECTRUM 

Calculations were performed for a s M = 1.2, 1.5, 5, 6 and 00, where M represents the 
bare heavy quark mass in the NRQCD action. Figure ^| shows the simulation energy of 
the ground state as a function of a s M for k = 0.24. The huge 0(1/M 3 ) effect at smaller 
M values is due to the vacuum expectation value of the cio term. This large correction 
to the unphysical simulation energy does not discredit the convergence of NRQCD, but 
special care must be taken to ensure that the large vacuum value is incorporated into the 
heavy quark propagation appropriately. In particular, previous work j|] has shown the linear 
approximation to be insufficient for the c\q term, which contains the vacuum expectation 
value, in computations of the S-wave spin splitting in the charm region. Fig. |2] explicitly 
shows the error introduced by placing the C\q term in Hb rather than Ha in Eq. (^). 
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l/[aM] 

FIG. 2. The simulation energy of a ground state heavy- light meson at rest. Results are displayed 
from terms up to 0(l/(a s M) k ), with k = 1, 2, 3. k is fixed at 0.24, and M is the bare heavy quark 
mass. Solid circles are k = 1, solid squares are = 2, open circles (squares) are k = 3 with the cio 
term in Hb(Ha) of Eq. (||). Statistical errors are smaller than the plotted symbols. 



It will also be noted from Fig. that the simulation energy is negative for a s M = 1.2 
at third order in NRQCD. Of course the absolute energy scale is unphysical in NRQCD 
due to omission of the large heavy quark mass term from the leading order action, and 
physical quantities (i.e. mass differences) are independent of the absolute energy scale. If 
the simulation energy were large and negative, it might signal a poor 1/M expansion and/or 
a problem for the discretization of heavy quark propagation, but a small negative result 
presents no problem. 

Table |IT| shows a calculation of the kinetic mass, and indicates that the bottom quark 
requires a s M « 5.5. The charm quark seems to want 1.2 < a s M < 1.5, although the data 
suggest that the 1/M expansion might not be converging in this region. It is interesting to 
notice that the vacuum expectation value does not affect the calculation of this observable 
significantly. 

A stronger statement about convergence comes from the splitting between the spin- 
singlet and spin-triplet S-wave states, since the uncertainties are smaller. Table [TV] suggests 
a nice convergence in the bottom region, but no guarantee of convergence for charm. It will 
be noted that an incorrect treatment of the vacuum expectation value (i.e. putting the c w 
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TABLE III. The kinetic energy of a 1 S'o heavy-light meson. and Hb are defined by Eq. 
@ and cio by Eq. (p"5|). Except for a s M = 6, the results in physical units are computed from 
0(1/M 2 ) data, using the lattice spacing from Eq. ( |30| ) to set the physical length scale. Only 
statistical uncertainties are shown. 



a s M 




a t E(p) 


- a t E(6) 




M kin [GeV] 




0(l/M) 


0(l/M z ) 


0(1/M 3 ) 










cio in F B 


cio in fi^ 




k=0.23 1.2 


0.0487(10) 


0.0515(10) 


0.0530(10) 


0.0582(39)* 


1.83(6) 


1.5 


0.0425(9) 


0.0444(9) 


0.0459(9) 


0.0519(25)* 


2.12(8) 


5.0 


0.0185(10) 


0.0188(10) 


0.0188(10) 


0.0239(26)* 


5.0(3) 


6.0 


0.0162(10) 








5.8(4) 


k=0.24 1.2 


0.0496(25) 


0.0523(24) 


0.0549(24) 


0.0566(27) 


1.80(10) 


1.5 


0.0430(24) 


0.0449(22) 


0.0469(21) 


0.0471(23) 


2.10(12) 


5.0 


0.0191(17) 


0.0194(16) 


0.0194(16) 


0.0194(16) 


4.9(4) 


6.0 


0.0180(22) 








5.2(7) 



*These computations use only 200 configurations. 



term into Hb) can actually lead to a small 0(1/M 3 ) contribution, but this is incidental. 

Tables [TLT| and [TV| correctly accommodate the vacuum expectation value of the cio term by 
placing it in Ha- In Ref. [|J, this method was found to give the same numerical results, within 
statistical uncertainties, for the S-wave kinetic mass and spin splitting as was obtained by 
computing the vacuum expectation value directly and subtracting it from the Hamiltonian. 
For the present calculation, a similar check was performed: the vacuum expectation value 
was computed from 400 of the gauge field configurations, and the S-wave kinetic mass and 
spin splitting were computed from 200 configurations using the Hamiltonian with the vacuum 
value explicitly removed. As expected, the results agree within statistics with Tables [TTT] and 
|TV] when the cio term is in Ha- 

Interpolating these data so that is the physical mass of a bottom meson, one finds a 
spin splitting which is only about 55% of the experimental value. This is typical of quenched 
lattice calculations (see for eg. Refs. 0,^,f2^|). Even an unquenched NRQCD calculation 
did not reproduce the experimental B* — B splitting [28|| so perhaps the tadpole-improved 
classical values, which were used for the coefficients q in the NRQCD Hamiltonian, account 
for the residual discrepancy. 

The data reported in Ref. [|J display a 1/M expansion for charmed mesons in lattice 
NRQCD in which 0(1/M 3 ) terms were as significant as 0(1/M 2 ) terms, so convergence of 
the expansion could not be assured. In that work, it was hoped that a more convergent 
expansion might be obtained via changes in the lattice method. In particular, replacement 
of the average plaquette tadpole factor by the mean link in Landau gauge was suggested 
to hold some promise. |||7| The present work has made this modification plus a number 
of others including: a more aggressively-improved light quark action, asymmetric lattices 
with temporal spacing reduced by a factor of two, Dirichlet temporal boundaries for light 
quarks rather than periodic ones, smeared meson sources rather than local sources, and a 
symmetric dependence on Hb in Eq. @. Despite these changes, convergence of the 1/M 
expansion remains uncomfirmed for S-wave charmed mesons. It is possible that a study of 
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TABLE IV. The 3 Si - 1 S mass splitting. H A and H B are defined by Eq. © and c w by Eq. 
(|15|). Except for a s M = 6, the results in physical units are computed from 0(1/M 2 ) data, using 
the lattice spacing from Eq. ( |30| ) to set the physical length scale. Only statistical uncertainties are 
shown. 



a s M 




a t M( d 5i) 


- atMC-So) 




M( d Si) - M^S'o) 




0(1/M) 


0(l/M z ) 


(9(1/M 3 ) 


[MeV] 








cio in 


cio in 




k=0.23 1.2 


0.0441(6) 


0.0527(7) 


0.0508(7) 


0.0566(25)* 


96.6(24) 


1.5 


0.0385(6) 


0.0447(6) 


0.0447(6) 


0.0452(19)* 


81.9(21) 


5.0 


0.0155(4) 


0.0165(4) 


0.0165(4) 


0.0159(8)* 


30.2(10) 


6.0 


0.0136(5) 








24.9(11) 


At=0.24 1.2 


0.0483(11) 


0.0589(12) 


0.0577(11) 


0.0677(15) 


107.9(32) 


1.5 


0.0418(9) 


0.0494(11) 


0.0502(10) 


0.0554(11) 


90.5(28) 


5.0 


0.0147(13) 


0.0157(13) 


0.0158(13) 


0.0158(13) 


28.8(25) 


6.0 


0.0136(7) 








24.9(14) 



*These computations use only 200 configurations. 



terms beyond 0(1/ M 3 ) would reveal that the series really is well-behaved, but this remains 
unexplored. 



V. P-WAVE SPECTRUM 



Each of the four P-wave operators from Table | leads to a visually-identifiable plateau; 
an example is shown in Fig. |3|. The method of maximum Q, discussed in Sec. |T|, can 
be used to define precise plateau boundaries and the resultant 3 P — 1 5'o splitting is shown 
in Table |V|. In contrast to the S-wave splitting discussed in the previous section, the 1/M 
corrections to the 3 Pq — 1 Sq splitting are not large relative to the statistical uncertainties, 
even in the charm region. The splittings given in Table [V] are consistent with the available 
experimental data, as will be discussed below. 

Tables and |V11| contain lattice results for splittings which involve the other P-wave 
mesons. In the charm region, the k = 0.23 data produce a 3 P2 meson which is heavier 
than the 3 Po meson, but at k — 0.24 the 3 p2 — 3 Po splitting is consistent with zero. The 
magnitudes of the splittings decrease in the bottom region, as expected. 

A comparison of Tables 0, [VI] and [VIl] plus the effective mass plots in Fig. |] provide some 
indication of the systematic uncertainty which arises from the choice of plateau region. In 
particular, it will be noted that the effective mass plots are monotonically decreasing near 
the source, so if a plateau region is chosen too near the source, it will produce a mass 
splitting which is too large. Thus one arrives at an upper bound for the 3 P2 — 3 Po splitting, 
as presented in Ref. ||24fl . 



The correlation functions constructed using 1 Pi and 3 Pi operators contain some combi- 
nation of the physical J = 1 mesons. Both calculations should lead to the same (lighter) 
physical mass at large Euclidean times if both operators have a substantial overlap with the 
less massive J = 1 state. In principle, the masses of the physical states can be obtained 
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TABLE V. The 3 P - mass splitting. H A and H B are defined by Eq. © and c 10 by 
Eq. (|T5|), Except for a s M = 6 and a s M = oo, the results in physical units are computed from 
0(1/ M 2 ) data, using the lattice spacing from Eq. ( |30| ) to set the physical length scale. Only 
statistical uncertainties are shown. 



a s M 






- a t M( l S Q ) 








Oi\ IM\ 


Oil /M 2 ) 


0(1/M 3 ) 


[MeVl 








cio in H B 


cio in H A 




k=0.23 1.2 


0.275(5) 


0.271(4) 


0.276(4) 


0.277(16)* 


497(13) 


1.5 


0.266(5) 


0.263(4) 


0.267(4) 


0.268(10)* 


482(13) 


5.0 


0.238(5) 


0.238(4) 


0.238(4) 


0.252(11)* 


436(12) 


6.0 


0.236(5) 








432(13) 


oo 


0.225(5) 








412(13) 


k=0.24 1.2 


0.311(6) 


0.311(6) 


0.315(6) 


0.321(7) 


570(16) 


1.5 


0.300(6) 


0.301(6) 


0.304(6) 


0.306(6) 


552(16) 


5.0 


0.256(8) 


0.256(7) 


0.256(7) 


0.256(7) 


469(16) 


6.0 


0.254(8) 








465(18) 




0.210(20) 








385(38) 



*These computations use only 200 configurations. 



TABLE VI. The 1 P 1 - ^o, 3 Pi - 1 S and 3 P 2 - mass splittings. Data at a s M = 6 use 



the Hamiltonian up to 0(1/M) and data at a s M < 6 use the Hamiltonian up to 0(1 /M 2 ). Only 
statistical uncertainties are shown. 





a s M 




a t M(X) - atM^So) 








X = l P x 


X = 3 Pi 


X = 3 P 2 


k=0.23 


1.2 


0.288(7) 


0.298(9) 


0.325(7) 




1.5 


0.279(7) 


0.285(8) 


0.311(7) 




5.0 


0.236(8) 


0.246(7) 


0.259(7) 




6.0 


0.240(7) 


0.243(8) 


0.253(7) 




oo 


0.181(26) 


0.226(10) 


0.214(11) 


k=0.24 


1.2 


0.295(13) 


0.309(20) 


0.327(18) 




1.5 


0.295(10) 


0.294(19) 


0.315(17) 




5.0 


0.232(15) 


0.250(14) 


0.270(13) 




6.0 


0.225(16) 


0.246(15) 


0.265(13) 




oo 


0.200(32) 


0.257(61) 


0.219(17) 
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FIG. 3. The simulation energy of a 3 Po heavy-light meson, for k = 0.23 and a s M = 5.0. The 
Hamiltonian contains all terms up to and including 0(l/(a s M) 2 ). The splitting between the 3 Po 
and 1 Sq is also shown. 



using the 1 P±/ 3 Pi operator basis by also calculating the mixing matrix elements, but the 
effect was too small to be observed in these data. In practice, for Euclidean times which 
can be used in our lattice simulation, no significant energy difference is observed between 
the 1 P\ and 3 Pi channels. 

Predictions for the physical mesons are displayed alongside experimental data in Table 
VIII] . The D s — D and B s — B mass differences depend only mildly on 0(1/M 2 ) and 
0(1/M 3 ) terms, and they are close to experiment. In contrast, the S-wave spin splittings 
are significantly smaller than experiment (as is typical of quenched lattice QCD |4],[23|,f25|j). 
Furthermore, the 0(1/M 2 ) and 0(1/M 3 ) corrections to D* — D are each about 20% of the 
leading order result. The results in Table [VIII| use the value of k s determined from uik in 
Eq. (|29|). Use of the other determination in Eq. ( p9t ) shifts the central values of the D* — D s 
and B* — B s mass splittings by 1 MeV or less, and shifts the splittings among P-waves (D** 
and B**) by 7 MeV or less. 

A detailed comparison of P-wave results would be somewhat premature, since the ex- 
perimental data are rather incomplete and often rely on theoretical models for input, while 
the lattice calculation is quenched and lacks a firm connection between the physical mesons 
and the J = 1 operators. Nevertheless, Table |VIII| shows a general consistency between the 
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TABLE VII. The x Pi - 3 P , 3 Pi - 3 P and 3 P 2 - 3 Po mass splittings. Data at a s M = 6 use 



the Hamiltonian up to 0(1/M) and data at a s M < 6 use the Hamiltonian up to 0(1 /M 2 ). Only 
statistical uncertainties are shown. 





a s M 




a t M(X) - a t M( 3 P ) 








X = 1 P\ 


X = 3 Pi 


X = 3 P 2 


k=0.23 


1.2 


0.017(7) 


0.027(9) 


0.054(7) 




1.5 


0.015(7) 


0.021(8) 


0.047(7) 




5.0 


-0.002(8) 


0.008(6) 


0.021(7) 




6.0 


0.004(7) 


0.007(7) 


0.017(7) 




oo 


-0.044(26) 


0.001(9) 


-0.011(11) 


k=0.24 


1.2 


-0.016(14) 


-0.002(20) 


0.016(19) 




1.5 


-0.005(10) 


-0.007(19) 


0.014(17) 




5.0 


-0.024(16) 


-0.006(13) 


0.014(14) 




6.0 


-0.029(17) 


-0.008(14) 


0.011(15) 




oo 


-0.010(32) 


0.047(21) 


0.009(27) 



experimental and computed P-wave masses. 

It is instructive to compare lattice P-wave masses to the predictions of models, such as 



quark models |14| , |20| , |21[| , a Bethe-Salpeter study [pL5|| , a chromodynamic potential model |L6 



a bag model |I7j and a Blankenbecler-Sugar approach ||18|| . Many of these present results for 
J = 1 directly in the l Pi/ 3 P\ basis which simplifies the comparison to lattice data, and the 
quark models may in general be more closely related to the quenched approximation than 
to experiment. 



The models of Refs. [|l4]- |il|[T8| predict the traditional hydrogen-like ordering of P- waves, 
where the 3 Pq is the lightest meson, the 3 P 2 is the heaviest, and the Pi states lie in between. 
The authors of Ref. fI7[] find, from lightest to heaviest, 3 P , A(3/2), 3 P 2 , Pi(l/2), where 



the arguments represent the angular momentum of the light degrees of freedom. Ref. f20 
predicts a dramatic inversion where the 3 Po is heavier than the 3 P 2 by 100(150) MeV for 
£)**(B**) mesons, but the Pi states are the absolute lightest and heaviest. Ref. pl[| claims 
very small splittings (tens of MeV) where the P(l/2) and P(3/2) doublets overlap to produce 
different orderings for the D, D s and B systems (the B s is ordered like the D). 

Fig. H shows the D\ — Dq and B\ — B^ splittings from lattice QCD and from the models 
just mentioned. Despite the range of model predictions, it is evident that some general 
consistency exists between our results and a number of the models. Notice in particular that 
our results are numerically distinct from the large inversion of Ref. [2D]. Fig. |BJ compares 



the 3 Po — 1 So splittings as obtained from lattice QCD and the models. 

Of special importance is the comparison with Ref. p5| , where the spectrum was also 
computed from quenched lattice NRQCD. The discrepancy between the two lattice calcu- 
lations is exemplified by Fig. |5|. While there are many differences in method between the 
two computations, it is difficult to identify a compelling reason for the disagreement. Fig. 
f| indicates that our data must satisfy M(P4) — M(Pq) < 100 MeV for any chosen plateau 
region, and this is not consistent with Ref. P5| . It is hoped that further lattice efforts will 
improve this situation. Unfortunately, a recent lattice NRQCD study of the heavy-light me- 
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TABLE VIII. The heavy-light spectrum compared to experiment. The hadron naming scheme 
of Ref. [^9| is followed |3(|. The charm and bottom masses are fixed to the experimental D s and 
B s masses (giving a s M c = 1.37(8) and a s M^ = 5.8(lg 3 )). The first uncertainty on lattice data 



combines the statistical error with the uncertainty of at from Eq. (30). The second uncertainty on 
lattice data corresponds to the errors in a s M c and a s M^. The value of k s determined from itik in 



Eq. (29) has been used. 



charmed meson masses [Me VI 



bottom meson masses [MeV] 







experiment /Kef. 


lattice 






_i_ /T"» j? 

experiment / Kef. 


lattice 






1969(1)/ H 


input 


B s 




5369(2)/ JH 


input 


D* s - 


D s 


144/ m 


0(1/M):79(2)(3) 
0(1/M Z ):94(2)(5) 
O(l/M 3 ):103(3)(6) 




B s 


47(3)/ @ 


26(1)(1 3 ) 


D* s0 - 


D s 




530(13)(5) 




B s 




451(15)(± 3 ) 


D' sl - 


D s 




1 Pi:531(16)(2) 


B' sl - 


B s 




1 P 1 :425(24)(i 2 ) 


D 8 i- 


D s 


566(1)/ H 


3 Pi:542(22)(8) 


B s i — 


B s 




3 P i: 449(22)(^) 




D s 


604(2)/ || 


585(19)(6) 


B s2 ~ 

(B? - 


B s 
B s ) 


484(15)/ @t 


478(21)(1« 2 ) 


D s - 


D 


99,104/ (2| 


O(l/M):105(2)(l) 
O(l/M 2 ):107(3)(l) 
0(1/M 3 ):112(4)(1) 


B s - 


B 


90(2)/ H 


92(3)C5) 


D* - 


D 


141,142/ [29| 


0(1/M):84(3)(3) 
O(l/M 2 ):101(3)(5) 
0(1/M 3 ):117(4)(6) 


B* - 


B 


46/ @ 


25(2)(1 3 ) 


D* - 


D 




579(15)(5) 


B*o~ 


B 




475(19)(±*) 


D[- 


D 


596(53)/ |ig 


1 Pi:548(20)(2) 


B[- 


B 


391(16)/ @* 


1 Pi:407(38)(l^) 


Di- 


D 


558(2)/ 


3 Pi:557(32)(8) 


Bt- 


B 


431(20)/ Q* 


3 P i: 453(33)(l? ) 


D* 2 - 


D 


595(2)/ H 


588(30) (6) 


B* 2 - 
(B** - 


B 

-B) 


489(8)/ H* 
460(13)/ @* 
418(9)/ ||f 


493(29)(^ 3 ) 



f Experimental signal is a sum over resonances with differing momenta J = 0, 1, 2. 
*Theoretical estimates for some of the mass splittings have been used as input. 
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son spectrum |26] has statistical uncertainties which are too large to resolve the discrepancy 
between our results and those of Ref. [[25 . 



VI. CONCLUSIONS 

The masses of S and P-wave heavy-light mesons have been calculated in the quenched 
approximation, using lattice NRQCD for the heavy quark and a highly-improved action for 
the light degrees of freedom. Calculations at first, second and third order in the heavy 
quark mass expansion were used as a test of convergence, and it was concluded that the 
0(1/M 2 ) and 0(1/M 3 ) contributions to the D* - D splitting are both near 20%. The 
D* — D s splitting receives a 20% correction at 0(1/M 2 ) and a 10% correction at 0(1/M 3 ). 
These results might represent a convergent 1/M expansion if terms beyond 0(1/M 3 ) are 
decreasing appropriately, but convergence cannot be ascertained from the present study. 
This same conclusion was reached in Ref. || by a computational method which differed 
from the present one in some details; most notably, Ref. || used the average plaquette 
tadpole definition whereas the present work uses the mean link in Landau gauge. 

No convergence problem is found for P-wave charmed masses, and the P-wave spectrum 
for both charmed and bottom mesons is predicted. The 3 P2 is heavier than the 3 Pq and P\ 
states in the D** system, with M(D* 2 ) - M(D* s0 ) = 55 ± 17 MeV. For the D**, B** and B** 
systems the P2 — Pq splittings are smaller and consistent with zero. These conclusions agree 
with a number of model calculations and are compatible with the available experimental 
data. 
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FIG. 4. Effective mass plots for the 3 P2 — 3 -Po mass splitting with k = 0.23 and a s M = 1.5, 5.0 
and oo. The plateau region and value is also shown. 
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FIG. 5. The D* 2 - DjJ and B\ - B$ splittings (solid symbols) and the D* 2 - D* and B* 2 - B* { 
splittings (open symbols) from lattice QCD and various model calculations. 
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FIG. 6. The Dq — D and Bq — B splittings (solid symbols) and the D* s0 - D s and B 
splittings (open symbols) from lattice QCD and various model calculations. 
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